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The paper investigates non- stationary heat conduction in one-dimensional models with substrate 
potential. In order to establish universal characteristic properties of the process, we explore three 
different models — Frenkel-Kontorova (FK), phi4+ (</> 4 +) and phi4- (4> 4 — ). Direct numeric sim- 
ulations reveal in all these models a crossover from oscillatory decay of short-wave perturbations 
of the temperature field to smooth diffusive decay of the long-wave perturbations. Such behavior 
is inconsistent with parabolic Fourier equation of the heat conduction and clearly demonstrates 
the necessity of hyperbolic models. The crossover wavelength decreases with increase of average 
temperature. The decay patterns of the temperature field almost do not depend on the amplitude 
of the perturbations, so the use of linear evolution equations for temperature field is justified. In 
all model investigated, the relaxation of thermal perturbations is exponential - contrary to linear 
chain, where it follows a power law. However, the most popular lowest-order hyperbolic general- 
ization of the Fourier law, known as Cattaneo-Vernotte (CV) or telegraph equation (TE) is not 
valid for description of the observed behavior of the models with on-site potential. In part of the 
models a characteristic relaxation times exhibit peculiar scaling with respect to the temperature 
perturbation wavelength. Quite surprisingly, such behavior is similar to that of well-known model 
with divergent heat conduction (Fermi-Pasta-Ulam chain) and rather different from the model with 
normal heat conduction (chain of rotators). Thus, the data on the non-stationary heat conduction 
in the systems with on-site potential do not fit commonly accepted concept of universality classes 
for heat conduction in one-dimensional models. 



I. INTRODUCTION 



Relationship between empiric equations of heat con- 
duction (Fourier law) and microstructure of solid di- 
electrics is known to be one of the oldest and most elusive 
unsolved problems in solid state physics [1, 2]. The clas- 
sic solution for the problem suggested by Peierls [3] was 
questioned after seminal numeric experiment of Fermi, 
Pasta and Ulam [4]; the latter has disproved common 
belief concerning fast thcrmalization and mixing in non- 
integrable systems with weak nonlincarity. Despite large 
amount of work done, necessary and sufficient conditions 
for microscopic model of solid to obey macroscopically 
the Fourier law with finite and size - independent heat 
conduction coefficient [4-16] are not known yet. Numer- 
ous anomalies in the heat transfer in microscopic models 
of dielectrics were revealed by means of direct numeric 
simulation over recent years, including qualitatively dif- 
ferent behavior of models of different types (with and 
without on-site potential) and dimensionality [1]. To 
date, it is believed that in one dimension the microscopic 
models with conserved momentum the heat conduction 
coefficient diverges in the thermodynamic limit (as the 
chain length N goes to infinity) as n ~ ./V 7 with 7 vary- 
ing in the interval 0.3-^0.4 [1]. The only known exception 
with convergent heat conduction coefficient is the chain 
of coupled rotators [8, 9]. As for the models without the 
moment conservation, their heat conduction is conver- 
gent [10, 14], with exception of integrable models [11]. In 
two dimensions, the divergence still exists [1], although is 
reported to be of logarithmic rather than power-like type. 
It is believed that in three dimensions the heat conduc- 



tion coefficient will converge, although some alternative 
data exist [15, 16]. 

Vast majority of results obtained to date in the field 
of microscopic foundations of the heat conduction dealt 
with stationary problem, with steady heat flux under 
constant thermal gradient. In all these cases the macro- 
scopic law to be checked was standard Fourier law of 
heat conduction. However, it is well-known that due to 
its parabolic character it implies infinite speed of the sig- 
nal propagation and thus should be modified if very large 
gradients or extremely small scales are involved [17-20]. 
On macroscopic level, numerous modifications were sug- 
gested to recover the hyperbolic character of the heat 
transport equation [19-22]. Perhaps, the most known 
is the lowest - order approximation known as Cattaneo- 
Vernotte (CV) or telegraphist (TE) law [23, 24]. It is 
written as 
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where k is standard heat conduction coefficient and r is 
characteristic relaxation time of the system, q is the vec- 
tor of heat flux and T is the temperature. The relaxation 
time is rather short for majority of materials (10 -12 s at 
room temperature), but can be of order 1 s, for instance 
in some biological tissues [25]. 

Modifications of the Fourier law bring about new ob- 
servable physical phenomena, such as temperature waves 
or the second sound [17-20]. Importance of the hy- 
perbolic heat conduction models for description of a 
nanoscale heat transfer has been recognized in recent ex- 
periments [26, 27]. 

Only few papers dealt with microscopic foundations of 
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non- Fourier heat conduction laws from the first principles 
or by its numeric investigation in the microscopic models 
[26, 28-37]. These works confirm that the non- Fourier 
effects may be very significant if large space gradients or 
fast changes of the temperature are involved. 

In recent paper [38] an attempt was made to relate 
the structure and parameters of microscopic models with 
conserved momentum to empiric description of the non- 
Fourier heat conduction. For this sake, two models be- 
longing to different universality classes with respect to 
the stationary heat conduction were studied: the /3-FPU 
chain and the chain of rotators. Oscillatory behavior of 
the decaying temperature disturbance, compatible only 
with hyperbolic macroscopic equation of the thermal 
transport, has been revealed in both models. The CV 
equation of the heat transport adequately describes the 
behavior of the chain of rotators, besides the region of 
crossover between the "hyperbolic" and the "parabolic" 
behavior. However, the /3-FPU model does not obey the 
CV equation. 

This paper continues the line of research started in 
[38] and investigates the relationship between structure, 
parameters and macroscopic description of the nonsta- 
tionary heat transfer in models with on-site potential. 
Such models are known to have the normal heat conduc- 
tion [9, 10, 14]. We are going to check (a) whether the 
effects of hyperbolicity can be observed in such models; 
(b) whether the thermal transport in these models can be 
described by linear equations and (c) whether the lowest 
- order CV equation is suitable for description of these 
phenomena. 



II. MODELS AND PROCEDURES 

The stationary heat transfer is characterized by sin- 
gle coefficient of the heat (or thermal) conductivity. In 
numeric experiments it is measured either by direct sim- 
ulation of the stationary flux under constant tempera- 
ture gradient, or from autocorrelation of the heat flux in 
isolated system via Green-Kubo formula [1]. The non- 
stationary conduction involves at least two parameters 
(CV equation) or even more in more advanced models. 
Moreover, one should not bind himself by a priori selec- 
tion of the empiric equation to fit the results. Then, it is 
desirable to find the characteristics of the process which 
can be measured directly from the numeric data and are 
not related to specific choice of the macroscopic empiric 
equation. 

One such choice may be observation of the temperature 
waves (the second sound), implemented, for instance, in 
papers [27, 35-37]. We adopt more general approach, 
based on relaxation profiles of different spatial modes of 
the temperature field. This approach is related to early 
numeric simulations on the non-stationary heat conduc- 
tion in argon crystals [30] and of thermal conductivity in 
super lattices [34]. 

In order to illustrate the approach, let us refer to ID 



version of the CV equation for temperature: 
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where a is the temperature conduction coefficient. 

To solve this equation, let us consider the problem 
of non-stationary heat conduction in a one-dimensional 
specimen with periodic boundary conditions T(L,t) — 
T(0,i), where T(x,t) is the temperature distribution, L 
is the length of the specimen, t > 0. If it is the case, one 
can expand the temperature distribution into Fourier se- 
ries: 



T(x,t) = cik(t) exp(27rifcx/L) 



(3) 
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with afc(i) = a,-k{t), since T(x,t) should be real, k is the 
modal number. 

Substituting (3) to (2), one obtains the equations for 
time evolution of the modal amplitudes: 
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Solutions of Eq. (4) are written as: 

ctk(t) = Ci k exp(Ait) + C 2 fe cxp(A 2 i), (5) 



where Ai >2 = (-1 ± y/l - m^War/I?) /2t ■, C lk and 
C2k are constants determined by the initial distribution. 

From (5) it immediately follows that for sufficiently 
short modes the temperature profile will relax in oscilla- 
tory manner: 



dk(t) ~ exp(— t/2r) exp(iwfct), 



(6) 



where u>k = \/16tt 2 k 2 ar/L 2 — 1/2t, k > L/Airy/ar. 

If the specimen is rather long [L 3> Att-Jot) then for 
small wavenumbers (acoustical modes): 
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(7) 



The first eigenvalue describes fast initial transient relax- 
ation, and the second one corresponds to stationary slow 
diffusion and, quite naturally, does not depend on r. This 
limit corresponds to Fourier heat conduction. 

So, we can conclude that there exists a critical length 
of the mode 



Iq = AiT\/aT. 



(8) 



which separates between two different types of the relax- 
ation: oscillatory and diffusive. 

The oscillatory behavior is naturally related to the hy- 
perbolicity of the system and is qualitatively inconsis- 
tent with the Fourier law. Here the critical scale is de- 
rived from the CV equation; however, it is clear that 
this scale, if it exists, can be revealed from numeric data 
without relying on any particular macroscopic descrip- 
tion. Moreover, this critical scale and its dependence on 
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the temperature and other parameters of the model help 
to understand how significant are the deviations from 
parabolicity in given model. In particular, if the size 
of the system is less than this critical mode length, the 
Fourier law cannot be used at all. 

Careful analysis of the relaxation profiles for different 
modes of the temperature perturbations also can give 
rather important insights. For instance, Equation (6) 
predicts that all oscillatory modes will decay with the 
same decrement related to unique relaxation time. One 
can verify this prediction by comparison with the nu- 
meric data and thus to evaluate the accuracy of the CV 
equation for particular model. 

In order to measure the critical size lo, the numeric ex- 
periment should be designed in order to simulate the re- 
laxation of thermal profile to its equilibrium value for dif- 
ferent spatial modes. For general one-dimensional model 
with on-site potential we simulate the chain of particles 
with unit masses with Hamiltonian 

N 

H = Y,[pI/ 2 + V(u n+ i - u n ) + U(u n )]. (9) 

n=l 

Here u n is the displacement of the n-th particle, p n = ii n , 
V(p) is the potential of the nearest-neighbor interaction 
and U(u) is the on-site potential (the minimum of V(p) 
corresponds to p = and minimum of U(u) is at it = 0). 
The global minimum of the potential energy corresponds 
to an unperturbed state {u n = 0}^ =1 . Boundary condi- 
tions are adopted to be periodic. 

In order to obtain the initial nonequilibrium tempera- 
ture distribution, all particles in the chain were initially 
connected to common Langevin thermostats. For this 
sake, the following system of equations was simulated: 

u n = V'{u n+ i - u n ) - V'(u n ~ u n -i) - U'{u n ) 

-iu n + £n, n = l...N, (10) 

where 7 is the damping coefficient and the white noise 
is normalized by the following conditions: 

(f n (t)) = 0, (£n(ti)U(h)) = 2iT n 5 nm 5{t l - t 2 ), (11) 

where T n is the prescribed temperature of the n-th par- 
ticle. 

In order to study the relaxation of different spatial 
modes of the initial temperature distribution, its profile 
has been prescribed as 

T n = T a + Acos[27r(n- l)/Z], (12) 

where Tq is the average temperature, A is the amplitude 
of the perturbation, and Z is the length of the mode 
(number of particles). The overall length of the chain N 
has to be multiple of Z in order to ensure the periodic 
boundary conditions. After the initial heating in accor- 
dance with (12), the Langevin thermostat was removed 
and relaxation of the isolated system to a stationary tem- 
perature profile is studied. The results were averaged 



over about 10 6 realizations of the initial profile {T n }^ =1 
in order to reduce the effect of fluctuations. 

In this paper, we analyze three models with linear 
nearest-neighbor interaction, which differ by the shape 
of the on-site potential. The Hamiltonian of all three 
models is written as: 

N 

H = X>n/ 2 + K+i - u n) 2 /2 + U(u n )]. (13) 

n=l 

Unit coefficient of the parabolic potential of the 
nearest-neighbor interaction does not effect the gener- 
ality. Real difference between the models appears due to 
different choices of the on-site potential U (u) . We treat 
three topologically different on-site potentials: sinusoidal 
potential (Frcnkcl-Kontorova model, FK) 

U(u) = 1-cosu; (14) 

(f> A - potential 

U{u) = 2u 2 (u-2ir) 2 /n 4 ; (15) 

</> 4 + potential 

U(u) = u 2 /2 + u 4 /4. (16) 

Substrate potentials (14), (15) and (16) differ topolog- 
ically. FK potential is periodic and bounded. Potential 
<j) 4 - is double- well. To simplify the comparison, a distance 
between the well minima and a height of the potential 
barrier are the same as for the FK model. Potential <fi + 
is single-well and unbounded. Our goal is to check how 
these differences affect the nonequilibrium heat transfer. 
From the viewpoint of equilibrium heat transfer, all these 
models have convergent thermal conduction coefficient in 
the thermodynamical limit [1]. 

III. RESULTS 

The first set of simulations reported here has been 
devoted to verification of transition from oscillatory to 
smooth relaxation profile of the temperature perturba- 
tion for different initial wavelength. Typical result of the 
simulation is presented in Fig. 1. The FK chains with the 
same length TV = 1024, average temperature To = 1 and 
amplitude of the perturbation A = 0.2 demonstrate qual- 
itatively different relaxation profiles for different modal 
wavelength Z . 

Results of the simulations presented in Fig. 2 demon- 
strate that all three investigated models with the on- 
site potential clearly exhibit similar transition from dif- 
fusive to oscillatory decay pattern as the characteris- 
tic wavelength decreases. It is possible to distinguish 
qualitatively three different relaxation patterns: diffu- 
sive (curves with k = 3), oscillatory (curves with k = 0, 
1) and crossover (curves with k = 2). In other terms, 
for all these models there exists the critical wavelength 




Figure 1: Relaxation of initial periodic thermal profile in the 
chain with FK potential, N = 1024, To = 1, A = 0.2 (a) 
Z = 32 and (b) Z — 128. Only part of the chain for n from 1 
to 128 is demonstrated. 

scale I* which separates the oscillatory from the diffusive 
decay. For the simulations presented here this character- 
istic scale is 

Z*~ 200-^ 300. (17) 

Thus, already at this step one can conclude that the 
"hyperbolic" behavior can be observed in all three mod- 
els. 

The wide range obtained in (17) is due to relatively 
small chain length; this simulation is sufficient to reveal 
the existence of transition between the decay patterns, 
but not sufficient for more or less exact estimation of the 
critical wavelength scale. Still, one can mention that sim- 
ilar values of the critical scale in Fig. 2 were obtained for 
different values of the average temperature of the system. 
Thus, one can conjecture that I* should be dependent on 
the temperature - similar observation was made also for 
the models with conserved momentum [38]. To verify 
this conjecture, we have performed more exact measure- 
ments of I* (for longer chains, up to 10000 particles, and 
with higher resolution) for all three models. The results 
are presented in Fig. 3. As the temperature decreases, 
the crossover length I* increases monotonously. 

IV. LINEARITY OF THE THERMAL 
RELAXATION PROFILES 

In order to check to what extent it might be possible 
to assume that the macroscopic equations of the nonsta- 
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Figure 2: (Color online) Evolution of the relaxation profile 
in the chain (length N = 800) for Z = 50, 100, 200 and 
400 (k = 0, 1, 2 and 3) with (a) FK potential (T = 4, 
A = 0.5); (b) 4 - potential (To = 4, A = 0.5); (c) (/> 4 + 
potential (T = 0.5, A = 0.05). Time dependence of the mode 
maximum T(l + Z/2) [red (gray) lines] and minimum T(l) 
[blue (black) lines] are depicted. Note that the time scales are 
different for each curve in order to fit them at one figure. 



tionary heat conduction for given models are linear, we 
have simulated the evolution of initial temperature distri- 
bution (12) with varying amplitude of the perturbation 
A, keeping all other parameters fixed. Then, the value 
AT n (t) = (T n (t)-T )/A is plotted versus time. If the re- 
sults coincide for different values of A, then it is possible 
to conclude that the relaxation of the thermal perturba- 
tions can be described by linear equation. In Fig. 4 the 
simulation results for the FK model arc demonstrated. 

For all regimes of the relaxation, one barely sees any 
difference for the relaxation profiles with different initial 
amplitudes; it happens despite the fact that the pertur- 
bation is by no means small - it achieves the half of the 
average temperature. Quite similar results were obtained 
also for (j) 4 - chain and 4 + chains and are omitted for the 
sake of brevity. So, one can conclude that the process of 
thermal relaxation in all these models may be described 
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Figure 3: Dependence of the critical wavelength I* on average 
temperature To of the chain with (f> 4 + substrate (curve 1), FK 
substrate (curve 2) and (j> 4 - substrate potential (curve 3). 



macroscopically with the help of linear equations with 
rather good accuracy. Similar conclusion was achieved 
also for the systems with conserved momentum, but with 
somewhat different method [38]. 

This conclusion does not seem to be trivial. The mod- 
els under consideration are essentially nonlinear and it 
is not completely clear why linear macroscopic equations 
should be so suitable for the description of the thermal re- 
laxation, especially for relatively large perturbation am- 
plitudes. 



V. THERMAL RELAXATION IN 
OSCILLATORY REGIME AND THE CV 
APPROXIMATION 

As it was mentioned before, the Cattaneo-Vernotte law 
predicts that the thermal relaxation will be exponential 
and that the relaxation time will be the same for all os- 
cillatory modes. Both these statements are not self - 
evident and should be verified numerically. 

For the sake of comparison, let us consider the same 
process of the relaxation of thermal waves in an infinite 
linear chain with the on-site potential. Such chain will be 
described by Equation (13) with the external potential 

U L (u) =wV/2. (18) 

Appropriate equations of motion will take the form 

u n + 2u n - u n -i - u n+ i + uj 2 u n = 0. (19) 

For the sake of simplicity, let us adopt that the initial 
temperature distribution is realized through attribution 
of initial velocities to each particle, with zero initial dis- 
placements. By virial theorem, such choice will not affect 
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Figure 4: (Color online) Dependence of rescaled temperature 
amplitude AT n (t)/A = (T n (t) - T )/A on rescaled time t/tj 
for FK chain: (a) Z = 50, tj = 0.5; (b) Z = 100, tj = 1; 
(c) Z = 200, tj = 2; (d) Z = 400, tj = 4. The length N = 
800, average temperature To = 4. Lines 1,3,5,7 [red (gray)] 
correspond to the maximum of the temperature in the chain 
n — 1, and Lines 2,4,6,7 [blue (black)] - to the minimum n = 
Z/2 + 1. Solid lines correspond to the perturbation amplitude 
A — 0.5, dashed lines - to A = 1, dots - to A — 2. The time 
scales are different for each curve in order to fit them at one 
figure. 



the long-time asymptotics. Namely, the initial conditions 
for Equation (19) are formulated as: 

u„(0) = 0, u n (0)=r] n , (r) n )=0, (20) 
(VnVm) = 2 Po + Acoa(2-Kn/Z)]S nm . 

where r\ n is random value with Gaussian distribution ac- 
cording to (20). The averaging is performed over the 
ensemble of the initial conditions. By standard trans- 
formations, it is easy to obtain the following solution of 
(19): 

OO 

«n(*)= X! VmG(n- m,t), (21) 

m— — oo 

G(p,t) = — I cos(t\/ oj 2 + 4 sin 2 \i — 2p^)dfx. 
Jo v 



The temperature distribution at arbitrary moment of 
time may be presented as: 

T n (t) = (u 2 n ) =T + 2Ag z {t) cos(27m/Z), (22) 

OO 

g z (t)= J2 G 2 ( P7 t)cos(2np/Z). 

P— — QO 

Time dependence of given thermal mode is completely 
governed by long-time asymptotics of function gz(t)- For 
given linear model with on-site potential it is difficult to 



6 



perform the exact integrations and summations in equa- 
tions (21), (22) (for the case without the on-site poten- 
tial these expressions can be reduced to compact exact 
expressions in terms of Bessel functions). Instead, we 
are going to derive the required asymptotics by approx- 
imate method, taking into account only generic features 
of the dispersion relation used in (21). For this sake, let 
us consider the general Green function similar to (21): 

G (p,t) = - f cos(</(a*) - 2pn)dfi, (23) 

where /(/i) is bounded at interval (0, n) and has a single 
inflection point /i* with positive first derivative, which 
satisfies the conditions 

/'V) =0, /V)>0, //e[0,7r]. (24) 

Then, for t large enough, the integral in (23) can be eval- 
uated by a stationary phase method. The condition of 
the stationary phase for specific values of p, t is presented 
as: 

Slfap, t) = tf(fx) - 2pn, O'(/i ) = tf'(no) -2p = 0. 

(25) 

For the sake of simplicity, we will consider only the 
case where (25) has unique solution. In the vicinity of 
this point, the phase is expanded as: 

A = f(no)t - 2pn Q + i/"0*o)t(p - Mo) 2 (26) 

+^rVoWM-Mo) 3 + .... 

6 

Standard evaluation of integral (23) for the case t ^s> p 
will, therefore, yield 

G (p,t) ~ r 1/2 , i->oo. (27) 

Estimation (27) is, however, wrong in the vicinity of 
the inflection point, which is defined by condition /iq = 
fx*. In this case, the term with the second derivative in 
expansion (26) will disappear and one will obtain: 

n = fOS)t - 2pp,* + V'VWm - n*) 3 + .... (28) 
6 

Quite obviously, evaluation of integral (23) with phase 
expansion (28) will yield: 

i r 

G (p,t) = - cos[fl(n,p,t)}dn 

7T Jo 

~ r 1 ' 3 cos[/(/i*)i - 2 P fi* + ip] + 0(r 2 ' 3 ), (29) 

where ip denotes constant phase shift. Expansion (29) is 
valid for 2p/t w /'(/x*), which, in turn, corresponds to 
maximum group velocity of linear waves in the system 
under consideration. Estimation (29) suggests, more ex- 
actly, that the expansion is valid for the interval of p 
defined as: 



Due to slower decrease rate with respect to time, the 
terms in expansion (22) in the interval of index defined 
by estimation (30) will be the most significant for long- 
time behavior. Then, the sum in (22) may be estimated 
as: 

p*+(p*) 1/3 2 
g k (t)~t- 2 / 3 £ cos(-^)cos 2 [/V*)t-2 w *+^ 

p—p* — (p* ) 1 / 3 

(31) 

For the most interesting case of relatively small 
wavenumbers the input of the sum in (31) will be of order 
of the square root of the number of participating sum- 
mands, i.e. of order (p*) 1 ^ ■ Consequently, with account 
of (30), one finally arrives to the following estimation: 

gz(t)-r^(p*) 1/6 -t-^. (32) 

As it is clear from the treatment, estimation (33) 
should be rather generic for linear chains. This power- 
law decay of the thermal mode is illustrated in Fig. 5, 

Therefore, we see that the thermal "relaxation" in lin- 
ear chains obeys power law rather than exponential. One 
can also mention that, quite as expected, the long-time 
asymptotics is not effected by the difference between the 
finite chain with periodic boundary conditions used for 
the simulations and the infinite chain used treated an- 
alytically. Also, it is not very surprising that the linear 
chains do not obey the CV law. However, one can expect 
that the nonlinear chains with on-site potential, known to 
have normal heat conduction [1] will obey it to some ex- 
tent. At least, the results for the models with conserved 
momentum [38] suggest such conjecture. 

The first statement to check is whether the decay of the 
thermal disturbances in the considered nonlinear models 
is indeed exponential. For this sake, we simulate the 
relaxation for different wavenumbers k for 4 + on-site 
potential (16). The results are presented in Fig. 6. This 
simulation clearly demonstrates the exponential charac- 
ter of decay for AT(t) = |Ti(i) - T \ + \T z/2+1 (t) - T \: 

AT(t) ~ exp(— t/t r ), when t -> oo, 

However in Fig. 6 we already can see that the char- 
acteristic relaxation time t r in the oscillatory regime de- 
pends on the perturbation wavelength, clearly at odds 
with the CV law (comp.(6)). To further clarify this point, 
we present in Fig. 7 the dependence of the characteristic 
relaxation time on the thermal perturbation wavelength 
for all three models with the on-site potentials studied 
above, as well as for the models with conserved momen- 
tum (FPU and chain of rotators) studied in [38]. 

As we can see from these results, three out of five mod- 
els (FPU, FK and ip 4 -) exhibit peculiar scaling of the re- 
laxation time in the oscillatory regime, which conforms 
to the approximate law 



\p-p*\~{p*)V\ p*=t/V)/2. (30) 



(33) 
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(a) 




Figure 5: (Color online) (a) Evolution of the relaxation profile 
in the linear chain (frequency cj = 1, length TV = 1024, aver- 
age temperature To = 2, amplitude A = 0.4) for Z — 32, 64, 
128 and 256 (j = 1, 2, 3 and 4; time scales ti = 1, ti = 1.5, 
t$ = 2 and £4 = 4 respectively). Time dependence of the 
mode maximum Ti [red (gray) lines] and minimum T z /2+i 
[blue (black) lines] are depicted. The time scales are different 
for each curve in order to fit them at one figure, (b) The 
same results are presented in a double logarithmic coordi- 
nates, AT(t) = |Ti(t)-T | + |T z/2+1 (f)-T |. Dashed fitting 
lines correspond to the power law AT ~ t -0,48 . 

with (3 w 0.45 for the FK and </- and (3 1.4 for the 
FPU model. As for ip 4 + model, the relaxation time also 
depends on the wavelength of the initial thermal profile, 
but it is difficult to speak about some peculiar scaling. 
Only in the chain of rotators the relaxation time in the 
oscillatory regime almost does not depend on the wave- 
length of the temperature perturbation. 



VI. CONCLUDING REMARKS 

On the basis of the simulations presented above, one 
can conclude that the one-dimensional models with topo- 
logically different on-site potentials exhibit qualitatively 
similar behavior with respect to nonstationary thermal 
conductivity. The thermal relaxation can be quite accu- 
rately described by linear equations even for relatively 
high perturbation amplitudes. All three models demon- 



(a) 




Figure 6: (Color online) (a) Evolution of the relaxation profile 
in the (f> 4 + chain (length TV = 1024, average temperature To = 

1, amplitude A = 0.2) for Z = 16, 32, 64, 128 and 256 (j = 1, 

2, 3, 4 and 5; time scale ti = 0.25, ti — ta = 0.5, ti = 1, 
and is = 2 respectively). Time dependence of the maximum 
Ti [red (gray) lines] and minimum T z / 2 +i [blue (black) lines] 
are depicted. The time scales are different for each curve 
in order to fit them at one figure, (b) Decay profile of the 
temperature difference AT{t) = |Ti (t) -T \ + \T z /2+i (i) - T 
in semilogarithmic coordinates. Dashed lines correspond to 
AT ~ exp(-t/t r ) with t r = 18.9, 22.2, 31.3, 100 and 370 for 
j = 1, 2, 3, 4 and 5 respectively. 



strate transition from oscillatory to diffusive relaxation 
regimes for growing wavelength of the initial thermal 
profile. Characteristic crossover wavelength rapidly de- 
creases with the temperature increase. 

Qualitatively, such behavior is compatible with the 
Cattaneo-Vernotte equation. However, quantitative 
study reveals that no unique relaxation time exists for 
different spatial harmonics of the initial temperature pro- 
file. For part (but not all) of the models, the relaxation 
times in the oscillatory regime approximately obey scal- 
ing law with respect to the wavelength of the initial ther- 
mal profile. 

All studied models are known to obey the Fourier law 
in the thermodynamical limit. The fact that the CV law 
is not suitable for description of the nonstationary heat 
transfer in these systems is rather surprising. Moreover, 
the scaling in FK and ip 4 - resembles that in the FPU 
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model, which has divergent heat conduction, although 
the scaling exponents are different. Still, at this stage one 
cannot claim that these exponents are in any sense uni- 
versal - a lot of further analysis is required. The chain of 
rotators turns out to be "exceptional" - it seems to obey 
the CV law both for the long and the short wavelengths 
of the temperature distribution. 



Figure 7: (Color online) Dependence of the characteristic re- 
laxation time t r on the period of initial thermal perturbation 
Z (double logarithmic scale): (a) the models with on-site po- 
tential: FK model (line 1, average temperature To = 2, per- 
turbation amplitude A = 0.5), (f> 4 — model (line 2, To = 2, 
A = 0.5), 4> i + model (line 3, T = 1, A = 0.2); (b) the mod- 
els with conserved momentum [38]: chain of rotators (line 4, 
T = 0.5, line 5, T = 0.3) FPU (line 6, T = 20). For all 
models, the solid part of the line corresponds to the oscilla- 
tory relaxation pattern and the dashed - to the crossover and 
monotonous relaxation. 



Scaling relationship (33) suggests that the macroscopic 
equations describing the nonstationary heat conduction 
in these models should include fractional derivatives. It 
is quite unusual that such "fractional' terms will appear 
as "hyperbolic" corrections to common Fourier law. 
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